# PBAR_ExamineBPData.py
#
# Examine borated-poly data
# 6/2013, John Kwong

from scipy import ndimage
import copy

# load files
baseDirectory = 'C:\Users\jkwong\Documents\Work\PBAR\data'
fileList = ['dx07', 'dx09', 'dx11', 'dw61', 'dw62','dw63','dw64','dw65','dw66','dw67', \
                     'dw68','dw69','dw74','dw75','dw76','dw77','dw78','dw79','dw80','dw81','dw82']

dat = []

for i in xrange(len(fileList)):
    fullFilename = baseDirectory + '\\' + fileList[i] + '.csv'
    print fullFilename
    
    dat.append(np.genfromtxt(fullFilename, delimiter = ','))

# dy98, read in differently
fullFilename = baseDirectory + '\\' + 'dy98Summed' + '.csv'
temp = np.genfromtxt(fullFilename, delimiter = ',')
temp2 = np.zeros_like(dat[0])
temp2[:,0] = dat[0][:,0]
temp2[:,1:] = temp.T
dat.append(temp2)
fileList.append('dy98')

# dz09, read in differently
fullFilename = baseDirectory + '\\' + 'dz09Summed' + '.csv'
temp = np.genfromtxt(fullFilename, delimiter = ',')
temp2 = np.zeros_like(dat[0])
temp2[:,0] = dat[0][:,0]
temp2[:,1:] = temp.T
dat.append(temp2)
fileList.append('dz09')

fileList = np.array(fileList)

# make groups
groups = dict()
groups['LessPolyReflective'] = np.array(('dw61', 'dw62','dw63','dw64','dw65','dw66','dw67','dw68','dw69'))
groups['MorePolyReflective'] = np.array(('dw74','dw75','dw76','dw77','dw78','dw79','dw80','dw81','dw82'))

keyNames = groups.keys()

groupsIndices = dict()
for i in xrange(len(keyNames)):
    temp = groups[keyNames[i]]
    temp2 = []
    for j in temp:
        temp2.append(find(fileList == j)[0])
    groupsIndices[keyNames[i]] = np.array(temp2)


# PLOT SPECTRA
# lining up spectra
# detector 23: normalize by counts and then corrected for gain

plotList = np.array((3, 6, 11, 12, 16, 21, 22))

detectorNumber = 35;
detectorNumber = 23;
detectorNumber = 9;


normalizeByCounts = 1
normalizeByPeakAmp = 0
normalizeManual = 0

manualNormalization= array([1, 1, 1, 1, 1, 0.0223, 0.0223])
manualNormalization = array([ 1.0, 0.74141671, 0.70887837, 0.78713622, 0.59239195, 0.45936951*0.0223*0.479, 0.44582909*0.0223*0.674])

filterSpectra = True
filterWidth = 2
correctGain = 1
correctGainManual = 0
#manualGainCorrection = array([ 1.0,  1.02535211,  1.00831025,  1.15,  1.15, 0.85, 0.85])
#manualGainCorrection = array([ 1.0, 1.0, 0.983, 1.119, 1.119, 0.824, 0.94])

countBound = [0.0035, 0.0050]
countBound = [0.0003, 0.0004]
#countBound = [0.0004, 0.0005]


plt.figure()
plt.grid()
j = 0
for i in plotList:
    
    x = dat[i][:,0][:]
    y = dat[i][:,detectorNumber][:]

    if filterSpectra:
        y = ndimage.filters.gaussian_filter(y, filterWidth, order = 0)    
    if normalizeByCounts:
        y = y/sum(y[1:-2])
    if normalizeByPeakAmp:
        cut = (x > 3.9) & (x < 4.1) 
        y = y / mean(y[cut])
    
    if normalizeManual:
        y = dat[i][:,detectorNumber][:] * manualNormalization[j]

    # adjust the gain
    if correctGain:
        if correctGainManual:
            x = x * manualGainCorrection[j]
        else:
            cut = (x > 1.0) & (x < 5.0) & (y > countBound[0]) & (y < countBound[1])
            x_mean = mean(x[cut])        
            x = x * 3.0 / x_mean
        
    plt.plot(x, y, label = fileList[i] + ", Det " + str(detectorNumber))
    j = j + 1

plt.xlim((0, 10))    
plt.yscale('log');
plt.legend()
plt.xlabel('Energy (MeV)');


# PLOT AVERAGE SPECTRA

keysList = ['LessPolyReflective', 'MorePolyReflective']

detectorNumber = 29;
plt.figure()
plt.grid()

for i in xrange(len(keysList)):
    
    y = np.zeros_like(dat[0][:,detectorNumber])
    
    for j in xrange(len(groupsIndices[keysList[i]])):
        y_temp = dat[groupsIndices[keysList[i]][j]][:,detectorNumber]
        y = y + y_temp/sum(y_temp)
    y = y / float(len(groupsIndices[keysList[i]]))
    x = dat[groupsIndices[keysList[i]][0]][:,0]
    
    plt.plot(x, y, label = keysList[i] + ", Det " + str(detectorNumber))

plt.xlim((0, 10))    
plt.yscale('log');
plt.legend()
plt.xlabel('Energy (MeV)');





if normalizeByCounts:
    plt.ylabel('Count (Normalized)');
else:    
    plt.ylabel('Count');


# PLOT SPECTRA, Difference

index1 = 1
index2 = 2
index2 = 0
detectorNumber = 21;

normalizeByCounts = True
normalizeByCounts = False


plt.figure()
plt.grid()

i = index2
x = dat[i][:,0][:]
y = dat[i][:,detectorNumber][:]

if normalizeByCounts:
    y = y/sum(y)
plt.plot(x, y, label = fileList[i] + ", Det " + str(detectorNumber))

i = index1
x = dat[i][:,0]
y = dat[i][:,detectorNumber]

if normalizeByCounts:
    y = y/sum(y)
plt.plot(x, y, label = fileList[i] + ", Det " + str(detectorNumber))

x = dat[index2][:,0]
y = dat[index2][:,detectorNumber] - dat[index1][:,detectorNumber]

if normalizeByCounts:
    y = y/sum(y)
plt.plot(x, y, label = fileList[index2] + " - " + fileList[index1] + ", Det " + str(detectorNumber))


plt.xlim((0, 10))    
plt.yscale('log');
plt.legend()
plt.xlabel('Energy (MeV)');


if normalizeByCounts:
    plt.ylabel('Count (Normalized)');
else:    
    plt.ylabel('Count');
